/* Copyright 2012 Tobias Marschall
 * 
 * This file is part of CLEVER.
 * 
 * CLEVER is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CLEVER is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CLEVER.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <iostream>
#include <fstream>
#include <math.h>
#include <boost/tokenizer.hpp>
#include <boost/lexical_cast.hpp>

#include "InsertLengthDistribution.h"

using namespace std;

InsertLengthDistribution::InsertLengthDistribution(const std::string& filename) {
	typedef boost::tokenizer<boost::char_separator<char> > tokenizer_t;
	boost::char_separator<char> separator(" \t");
	ifstream f(filename.c_str());
	if (f.fail()) {
		ostringstream oss;
		oss << "Could not open file \"" << filename << "\"";
		throw std::runtime_error(oss.str());
	}
	int last = -1;
	double psum = 0.0;
	string line;
	while (getline(f, line)) {
		tokenizer_t tokenizer(line,separator);
		vector<string> tokens(tokenizer.begin(), tokenizer.end());
		if (tokens.size()!=3) {
			throw std::runtime_error("Error parsing insert length distribution.");
		}
		try {
			int start = boost::lexical_cast<int>(tokens[0]);
			int end = boost::lexical_cast<int>(tokens[1]);
			double p = boost::lexical_cast<double>(tokens[2]);
			assert((last==-1) || (start==last+1));
			assert(start <= end);
			if (last==-1) min = start;
			max = end;
			last = end;
			psum += p;
			values.push_back(make_pair(start,p/(end-start+1)));
		} catch(boost::bad_lexical_cast &){
			throw std::runtime_error("Error parsing alignment pair.");
		}
	}
	if (fabs(1.0 - psum) > 1e-6) {
		throw std::runtime_error("Probabilities in insert size distribution do not sum to one");
	}
}

InsertLengthDistribution::~InsertLengthDistribution() {
}

double InsertLengthDistribution::probability(int insert_length) const {
	if (insert_length < min) return 0.0;
	if (insert_length > max) return 0.0;
	typedef vector<value_t>::const_iterator iter_t;
	iter_t cit = upper_bound(values.begin(), values.end(), make_pair(insert_length,2.0));
	--cit;
	return cit->second;
}

std::ostream& operator<<(std::ostream& os, const InsertLengthDistribution& ild) {
	typedef vector<InsertLengthDistribution::value_t>::const_iterator iter_t;
	for (iter_t cit=ild.values.begin(); cit!=ild.values.end(); ++cit) {
		os << cit->first << " " << cit->second << endl;
	}
	return os;
}
